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High-order  finite  element  methods  (also  known  as  spectral/hp  element  methods)  using  either  the  continuous  Galerkin  or  discontinuous 
Galerkin  formulation  have  reached  a  level  of  sophistication  such  that  they  are  now  commonly  applied  to  a  diverse  set  of  real-life 
engineering  problems.  Visualization  of  computed  results  is  often  used  as  a  means  of  understanding  and  evaluating  the  numerical 
approximation  of  the  mathematical  model,  and  it  provides  a  means  of  “closing  the  loop”  -  that  is,  of  critically  evaluating  the  computational 
results  for  refinement  of  the  model  and/or  numerics  or  for  interpretation  of  the  physical  world.  Visualizations  of  high-order  finite  element 
results  which  do  not  respect  the  a  priori  knowledge  of  how  the  data  were  produced  and  which  do  not  provide  a  quantification  of  the  visual 
error  produced  undermine  the  scientific  process  just  described.  The  goals  of  this  effort  are  to  define,  investigate,  and  address  the  technical 
obstacles  inherent  in  visualization  of  data  derived  from  high-order  numerical  methods  and  to  develop  algorithms  and  software  solutions  that 
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Abstract:  High-order  finite  element  methods  (also  known  as  spectral/hp  element 
methods)  using  either  the  continuous  Galerkin  or  discontinuous  Galerkin  formulation 
have  reached  a  level  of  sophistication  such  that  they  are  now  commonly  applied  to  a 
diverse  set  of  real-life  engineering  problems.  Visualization  of  computed  results  is  often 
used  as  a  means  of  understanding  and  evaluating  the  numerical  approximation  of  the 
mathematical  model,  and  it  provides  a  means  of  “closing  the  loop”  -  that  is,  of  critically 
evaluating  the  computational  results  for  refinement  of  the  model  and/or  numerics  or  for 
interpretation  of  the  physical  world.  Visualizations  of  high-order  finite  element  results 
which  do  not  respect  the  a  priori  knowledge  of  how  the  data  were  produced  and  which  do 
not  provide  a  quantification  of  the  visual  error  produced  undermine  the  scientific  process 
just  described.  The  goals  of  this  effort  are  to  define,  investigate,  and  address  the 
technical  obstacles  inherent  in  visualization  of  data  derived  from  high-order  numerical 
methods  and  to  develop  algorithms  and  software  solutions  that  can  be  employed  by  the 
high-order  simulation  community. 


Statement  of  Problem  Studies 

The  goals  of  this  effort  are  to  define,  investigate,  and  address  the  technical  obstacles 
inherent  in  visualization  of  data  derived  from  high-order  numerical  methods  and  to 
develop  algorithms  and  software  solutions  that  can  be  employed  by  the  high-order 
simulation  community. 


Summary  of  Results 

In  this  section,  we  first  present  the  motivating  work  (published  in  [Rl]  by  one  of  the 
investigators)  for  our  previous  ARO  grant  and  then  provide  a  summary  of  results  as  a 
consequence  of  ARO  funding.  To  summarize  -  six  peer-reviewed  journal  articles  have 
been  published  or  accepted  for  publication:  four  articles  targeting  the  visualization 
community  [J 1 ,  J2,  J4,  J5]  and  two  targeting  the  computational  mathematics  community 
[J3,J6]. 


Isosurface  Visualization 


Our  initial  work  on  high-order  finite  element  visualization  was  motivated  by  the  work  of 
Nelson  and  Kirby  [Rl],  in  which  they  presented  an  algorithm  for  ray-casting  high-order, 
spectral//?/)  elements.  Their  method  uses  a  world-space  approximation  of  the 
composition  of  the  coordinate  transformation  and  the  reference  space  basis  functions.  It 
assumes  multi-linear  mappings  (linear  element  boundaries  in  world  space),  and  includes 
a  quantification  of  the  approximation  and  root-finding  error.  They  show  that  the  image- 
space  method  compares  favorably  with  marching  cubes  in  compute  time  when  the 
tolerances  on  surface  position  are  sufficiently  high.  Figure  1  provides  an  example  of  the 
type  of  visualizations  produced  by  their  work.  The  marching  cubes  image  (left)  was 
generated  by  sampling  the  finite  element  volume  on  a  rectilinear  grid  of  spacing  $h$, 
using  a  marching  cubes  algorithm  to  provide  a  tessellated  isosurface,  and  rendering  the 
triangular  isosurface  using  ray-casting  (since  the  marching  cubes  result  is  a  triangular 
mesh,  the  ray-casting  can  be  done  exactly  as  done  in  [R2]).  For  the  marching  cubes 
image  presented,  a  grid  spacing  of  h=0.015  (yielding  4,705,274  voxels)  was  used.  For 
the  high-order  ray-traced  image  (right),  mapping  inversion  error  of  10'  and  11  order 
projected  polynomials  were  used.  These  parameters  were  chosen  such  that  the 
spectral//?/?  element  evaluation  time  and  rendering  time  was  nearly  identical  to  generate 
the  two  images.  The  root-mean-square  error  for  the  marching  cubes  image  is  0.0158;  the 
root-mean-square  error  for  the  ray-traced  image  is  3.5e-l  1.  The  images  look  very  similar, 
however  the  root-mean-square  error  difference  between  the  images  is  significant.  We 
should  also  point  out  that  the  file  size  for  the  marching  cubes  representation  is  over  an 
order  of  magnitude  larger  than  the  high-order  representation. 


Figure  1:  Marching  cubes  image  with  h=0.015  corresponding  to  4,705,274  voxels  (left) 
and  ray-traced  solution  using  11th  order  projected  polynomials  (right)  for  isosurface  of 
pressure  at  C  =  0.0  chosen  such  that  the  spectral//?/?  element  data  evaluation  and 
rendering  time  is  nearly  identical  (on  the  order  of  200  seconds).  The  root-mean- square 
error  for  the  marching  cubes  image  is  0.0158;  the  root-mean-square  error  for  the  ray- 
traced  image  is  3.5e-l  1. 


Although  the  ray-casting  methodology  provided  a  “pixel  exact”  visualization  of  the 
isosurface,  it  did  so  in  what  is  referred  to  as  “image  space”.  This  implies  that  even  after  a 
researcher  found  the  isosurface  of  interest  which  they  wanted  to  examine,  each  rotation, 
translation  or  zoom  into  the  image  required  approximately  the  same  amount  of  rendering 
time  as  each  pixel's  color  has  to  be  recomputed. 

The  classic  way  to  attempt  to  solve  this  issue  is  to  render  things  in  “object  space”  -  that 
is,  to  generate  objects  (triangles,  for  instance)  on  the  isosurface  so  that  once  an  isosurface 
is  found  and  an  object  is  created,  its  rendering  can  be  done  quickly.  In  [Jl]  we  proposed 
visualizing  isosurfaces  in  high-order  finite  element  datasets  with  a  particle  system  as  a 
means  of  solving  this  problem.  We  presented  a  framework  that  allows  particles  to 
sample  an  isosurface  in  reference  space,  avoiding  the  costly  inverse  mapping  of  positions 
from  world  space  when  evaluating  the  basis  functions.  The  distribution  of  particles  across 
the  reference  space  isosurface  is  controlled  by  geometric  information  from  the  world 
space  isosurface,  such  as  the  surface  gradient  and  curvature.  The  resulting  particle 
distributions  can  be  distributed  evenly  or  adapted  to  accommodate  world-space  surface 
features.  This  provides  compact,  efficient,  and  accurate  isosurface  representations  of 
these  challenging  data  sets.  In  Figure  2  we  present  a  visualization  of  an  isosurface  of 
pressure  within  an  incompressible  fluid  flow  field. 


Figure  2:  An  isosurface  of  a  finite  element  fluid  simulation  pressure  field  sampled  with  a 
particle  system.  The  color  indicates  the  relative  direction  of  the  surface  normal  at  the 
particle  (blue  indicates  outward  and  red  indicates  inward).  } 

When  one  employs  objects  to  mark  or  denote  an  isosurface,  one  faces  the  challenge  of 
knowing  how  many  objects  to  use  and  how  densely  to  pack  them.  A  sparse  packing  of 
the  objects  can  miss  critical  features  of  the  isosurface.  A  dense  packing  can  be  very 
inefficient  (especially  when  the  density  is  much  higher  than  is  needed).  In  [J2],  we 
describe  a  method  for  constructing  isosurface  triangulations  of  sampled,  volumetric, 
three-dimensional  scalar  fields  that  attempts  to  tackle  this  sampling  density  problem. 
The  resulting  meshes  consist  of  triangles  that  are  of  consistently  high  quality,  making 
them  well  suited  for  accurate  interpolation  of  scalar  and  vector-valued  quantities,  as 
required  for  numerous  applications  in  visualization  and  numerical  simulation.  The 
proposed  method  does  not  rely  on  a  local  construction  or  adjustment  of  triangles  as  is 
done,  for  instance,  in  advancing  wavefront  or  adaptive  refinement  methods.  Instead,  a 
system  of  dynamic  particles  optimally  samples  an  implicit  function  such  that  the 


particles'  relative  positions  can  produce  a  topologically  correct  Delaunay  triangulation. 
Thus,  the  proposed  method  relies  on  a  global  placement  of  triangle  vertices.  The  main 
contributions  of  this  work  was  the  integration  of  dynamic  particles  systems  with  surface 
sampling  theory  and  PDE-based  methods  for  controlling  the  local  variability  of  particle 
densities,  as  well  as  detailing  a  practical  method  that  accommodates  Delaunay  sampling 
requirements  to  generate  sparse  sets  of  points  for  the  production  of  high-quality 
tessellations.  In  [J5]  we  extended  this  work  to  handle  surfaces  that  come  as  a 
consequence  of  multi-material  interfaces. 

•  Streamline  Integration 

A  quick  search  of  both  the  visualization  and  the  application  domain  literature 
demonstrates  that  streamlines  are  a  popular  visualization  tool,  second  only  to  isosurfaces. 
The  bias  toward  using  streamlines  is  in  part  explained  by  studies  that  show  streamlines  to 
be  effective  visual  representations  for  elucidating  the  salient  features  of  the  vector  fields 
[R3].  Furthermore,  streamlines  as  a  visual  representation  are  appealing  because  they  are 
applicable  for  both  two-dimensional  and  three-dimensional  fields  [R4].  It  was  for  this 
reason  that  we  invested  time  considering  how  streamlining  would  be  impacted  by  high- 
order  finite  element  data. 

Streamline  integration  is  often  accomplished  through  the  application  of  ordinary 
differential  equation  (ODE)  integrators  such  as  predictor-corrector  or  Runge-Kutta 
schemes.  The  foundation  for  the  development  of  these  schemes  is  the  use  of  Taylor 
series  for  building  numerical  approximations  of  the  solution  of  the  ODE  of  interest. 

Taylor  series  can  be  further  used  to  elucidate  the  error  characteristics  of  the  derived 
scheme.  All  schemes  employed  for  streamline  integration  that  are  built  using  such  an 
approach  exhibit  error  characteristics  which  are  predicated  on  the  smoothness  of  the  field 
through  which  the  streamline  is  being  integrated. 

Low-order  and  high-order  finite  volume  and  finite  element  fields  are  among  the  most 
common  types  of  fluid  flow  simulation  datasets  available.  Streamlining  is  commonly 
applied  to  these  datasets.  The  property  of  these  fields  which  challenges  classic 
streamline  integration  using  Taylor  series  based  approximations  is  that  finite  volume 
fields  are  piecewise  discontinuous  and  finite  element  fields  are  only  C°  continuous. 
Hence  one  of  the  limiting  factors  of  streamline  accuracy  and  integration  efficiency  is  the 
lack  of  smoothness  at  the  inter-element  level  of  finite  volume  and  finite  element  data. 

Adaptive  error  control  techniques  are  often  used  to  ameliorate  the  challenge  posed  by 
inter-element  discontinuities.  To  paraphrase  a  classic  work  on  the  subject  of  solving 
ODEs  with  discontinuities  [R5],  one  must  (1)  detect,  (2)  determine  the  order,  size  and 
location  of,  and  (3)  judiciously  “pass  over”  discontinuities  for  effective  error  control. 
Such  an  approach  has  been  effectively  employed  within  the  visualization  community  for 
overcoming  the  challenges  posed  by  discontinuous  data  at  the  cost  of  increased  number 
of  evaluations  of  the  field  data.  The  number  of  evaluations  of  the  field  increases 
drastically  with  every  discontinuity  that  is  encountered  [R5].  Thus  if  one  requires  a 
particular  error  tolerance  and  employs  such  methods  for  error  control  when  integrating  a 


streamline  through  a  finite  volume  or  finite  element  dataset,  a  large  amount  of  the 
computational  work  involved  is  due  to  handling  inter-element  discontinuities  and  not  the 
intra-element  integration.  We  demonstrate  this  in  Figure  3  where  one  can  see  that  the 
number  of  streamline  sampling  steps  goes  up  drastically  each  time  a  streamline  attempts 
to  traverse  over  an  element  boundary. 


Figure  3:  The  center  graph  shows  a  streamline  on  an  $L_2$  projected  field  integrated 
using  RK-4/5.  The  left  graph  shows  the  streamline  between  t=0  and  t=0.3  and  the 
cumulative  number  of  RK-4/5  steps  (including  rejects)  required  for  integration.  Vertical 
lines  on  this  graph  represent  multiple  rejected  steps  occurring  when  the  streamline 
crosses  element  boundaries.  The  right  graph  shows  the  cumulative  number  of  RK-4/5 
steps  required  for  integration  to  t=2.0. 

As  the  root  of  the  difficulties  is  the  discontinuous  nature  of  the  data,  one  could  speculate 
that  if  one  were  to  filter  the  data  in  such  a  way  that  it  was  no  longer  discontinuous, 
streamline  integration  could  then  be  made  more  efficient.  The  caveat  that  arises  when 
one  is  interested  in  simulation  and  visualization  error  control  is  how  does  one  select  a 
filter  that  does  not  destroy  the  formal  accuracy  of  the  simulation  data  through  which  the 
streamlines  are  to  be  integrated?  Recent  mathematical  advances  [R6,  R7]  have  shown 
that  such  filters  can  be  constructed  for  high-order  finite  element  and  discontinuous 
Galerkin  (high-order  finite  volume)  data  on  uniform  quadrilateral  and  hexahedral  meshes. 
These  filters  are  such  that  they  have  the  provable  quality  that  they  increase  the  level  of 
smoothness  of  the  field  without  destroying  the  accuracy  in  the  case  that  the  “true 
solution”  that  the  simulation  is  approximating  is  smooth.  In  fact,  in  many  cases,  these 
filters  can  increase  the  accuracy  of  the  solution. 

As  part  of  our  work,  we  investigated  the  use  of  such  filters  applied  to  discontinuous  data 
prior  to  streamline  integration,  and  found  that  they  can  drastically  improve  the 
computational  efficiency  of  the  integration  process.  We  currently  have  two  published 
papers  on  this  topic  [J3,  J4]  (one  presenting  this  work  to  the  visualization  community, 
and  one  paper  presenting  new  computational  mathematics  work  which  came  as  a 
consequence  of  this  study).  We  also  have  an  accepted  paper  in  which  we  have  adapted 
this  idea  to  be  more  computationally  efficient  [J6] .  We  proposed  a  new  technique  that 
uses  a  one-dimensional  convolution  kernel  to  introduce  continuity  between  elements,  and 
increase  smoothness  while  not  introducing  additional  error  in  the  solution.  Furthermore, 
this  one-dimensional  implementation  is  the  same  regardless  of  the  dimension  of  the 


simulation  data.  This  in  turn  will  aid  in  accomplishing  the  goals  of  visualization  of  data 
over  more  complex  geometries  while  still  improving  the  smoothness  of  the  field  and  not 
compromising  the  accuracy  of  the  data. 
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